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ABSTRACT 

A  study  of  the  first  and  second-order,  dynamic,  ideal 
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are  obtained  with  a  FORTRAN  IV  computer  program  based  on  the 
finite  element  method  using  isoparametric  elements. 
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circular and  bulb  hull  form  including  finite  depth  effects 
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hull  form  in  heave  also  including  finite  depth  effects  are 
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I.   INTRODUCTION 

Naval  architects  and  marine  engineers  have  long  been 
greatly  interested  in  the  study  of  the  dynamic  fluid  response 
to  the  forced  oscillations  of  a  two-dimensional  floating  body 
The  applications  of  the  studies  include  studying  the  rigid- 
body  response  of  the  ship  due  to  an  incoming  wave  train,  and 
determining  the  natural  frequency  and  vibration  mode  shapes 
of  a  ship's  structure  in  connection  with  propulsion  systems 
design  [3^]. 

Marine  engineers,  in  modern  times,  are  concerned  with, 
among  other  things,  the  natural  frequencies  and  stability 
characteristics  of  floating  platforms  near  a  coastline. 
Also,  they  are  concerned  with  the  related  problem  of  finding 
the  dynamic  loads  on  submerged,  moored,  or  fixed  bodies  [8], 

The  modern  techniques  for  solving  these  types  of  problems 
involve  seeking  the  solution  of  a  steady-state,  periodic, 
potential  flow  problem  with  a  free  surface  in  a  fluid  of 
either  infinite  or  finite  depth.   This  problem  is  made 
tractable  by  linearizing  the  free  surface  boundary  condition. 

A.   HISTORICAL  BACKGROUND 

The  historical  development  of  potential  flow  problems 
involving  a  free  surface  with  wave  motion,  and  the  related 
problem  of  the  dynamic  response  of  fluids  due  to  the  motion 
of  rigid  bodies  on  and  under  a  free  surface  is  lengthy;  but 
of  interest.   Several  authors  have  reviewed  this  history  in 
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some  detail  [3,17,24,31].   The  history  is  repeated  in  the 
following  in  the  interest  of  providing  a  foundation  upon 
which  this  author  endeavors  to  add  additional  information  in 
an  attempt  to  provide  a  more  complete  history. 
1 .   First  Order  Solutions 

a.   General  Potential  Theory 

The  foundation  of  the  modern  day  theory  of  fluid 
motions  begins,  of  course,  with  the  classical  works  of  Newton, 
Lagrange,  Euler,  Stokes,  and  others. 

However,  in  1879,  Lamb  [15]  provided  one  of  the 
first  complete  studies  of  the  theory  of  wave  motions;  which 
he  updated  five  times  until  1932.  It  appears  that  the  next 
major  contribution  was  made  by  Havelock  [10]  in  1929;  wherein, 
he  studied  the  solution  of  the  wave  forms  created  by  a  wave 
generator. 

The  modern  beginnings  of  the  study  of  the  inter- 
action of  a  fluid  and  a  harmonically  oscillating  rigid  body 
at  the  free  surface  start  with  the  work  of  Lewis  [19]  in  1929. 
Lewis  studied  the  determination  of  added  mass  with  a  boundary 
condition  requiring  zero  pressure  at  the  mean  position  of  the 
free  surface,  thus  eliminating  the  effects  of  frequency  of 
body  oscillation  and  of  damping.   Lewis  considered  the  high 
frequency  case  of  the  more  general  problem. 

The  first  work  which  included  the  effects  of 
frequency  and  damping  was  done  by  Ursell  [38]  in  19^9.   This 
work  solved  the  problem  of  a  right-circular,  semi-immersed 
cylinder  heaving  on  the  free  surface  of  a  fluid  infinite  in 
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extent.   Ursell  presented  added  mass  and  damping  coefficients 
as  a  function  of  frequency.   His  method  of  solution  involved 
placing  an  "infinite"  number  of  singularities  of  all  orders 
at  the  intersection  of  the  free  surface  at  rest  and  the 
cylinder  centerline .   The  strengths  of  the  sources  were 
determined  from  the  kinematic  boundary  condition  applied  on 
the  immersed  surface  of  the  cylinder.   Grim  [9]  arrived  at  a 
similar  solution  in  1953. 

A  period  of  about  ten  years  passed  before  any 
additional  work  was  reported  on  the  subject.   There  is  a 
rather  obvious  reason  for  this  time  lag.   Ursell' s  solution 
involved  laborious  numerical  calculations  which  had  to  be 
done  by  hand.   The  advent  of  the  digital  computer  made  it 
feasible  to  pursue  further  solutions. 

In  I960,  Tasai  [36]  extended  the  work  of  Ursell 
to  include  roll  and  sway  motions  of  Lewis  forms.   At  the  same 
time,  Porter  [30]  compared  linearized  pressure  calculations 
with  experimental  results,  and  extended  the  range  of  solutions 
to  include  other  hull  forms. 

It  appears  that  Yu  and  Ursell  [42]  provided  the 
first  theoretical  study  of  the  effect  of  finite  depth  on 
two-dimensional  solutions  in  1961.   This  was  followed  by 
additional  work  by  Kim  [14]  in  1969,  who  also  included  finite 
depth.   In  the  first  work,  added  mass  coefficients  and  the 
properties  of  the  resultant  linear  wave  were  reported.   The 
second  work,  by  Kim,  gave  added  mass  and  damping  coefficients 
as  a  function  of  frequency  and  depth.   Work  by  Paulling  and 
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Richardson  [28]  and  Paulling  and  Porter  [29]  provided  some 
experimental  verification  of  the  theory  in  1962.   Vugts  [40], 
at  the  Netherlands  Research  Center  in  1968,  reported  additional 
extensive,  experimental  research  on  seven  different  ship  hull 
forms . 

b.   Finite  Element  Solutions 

Zienkiewicz  [46],  in  1965,  provided  the  first 
successful  application  of  the  finite  element  method  to  Poisson's 
equation.   Although  the  finite  element  method  was  originally 
developed  to  solve  problems  in  the  theory  of  elastic  mediums, 
Zienkiewicz ' s  work  "opened  the  door"  to  possible  solutions  to 
boundary  value  problems  of  many  different  types,  including  the 
solution  of  fluid  flow  problems;  this  follows  since  Laplace's 
equation  is  contained  in  Poisson's  equation. 

In  the  same  year  (1965) 9  Zienkiewicz,  Irons,  and 
Nath  [44]  first  used  the  finite  element  method  to  find  added 
nfass .   Just  as  in  the  case  of  Lewis  and  others  before,  surface 
waves  were  excluded  from  this  solution,  as  well  as  in  subse- 
quent works  by  Rftren  [32],  Holand  [12],  Matsuura  and  Kawakami 
[22],  and  Matsumoto  [21]. 

Zienkiewicz  and  Newton  [43]  presented  the  first 
general  theoretical  treatment  of  the  fluid-structure  inter- 
action problem  using  the  finite  element  theory  in  1969.   Their 
work  provided  a  means  of  including  surface  waves  and  deter- 
mination of  energy  loss  due  to  wave  propagation.   The  latter 
was  accomplished  by  providing  a  radiation  boundary  condition. 
This  radiation  boundary  condition  was  essential  to  the 
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successful  modeling  of  infinite  regions  under  steady-state 
conditions.  In  addition,  compressibility  effects  could  be 
included  in  the  analysis. 

Chenault  [4]3  in  1970,  appears  to  be  the  first  to 
apply  the  finite  element  analysis  developed  by  Zienkiewicz 
and  Newton  to  find  added  mass  and  damping  coefficients  as  a 
function  of  frequency  and  hull  form.   The  fluid  regions  he 
studied  were  chosen  to  simulate  infinite  depth.   Some  of  his 
results  are  published  in  Ref.  [23].   Bai  [3]  also  included 
surface  waves  in  his  1972  analysis  of  a  wide  range  of  problems. 
His  work  demonstrated  the  flexibility  of  the  method  by  studying, 
in  addition  to  the  problems  of  Chenault,  the  additions  of  sway 
and  roll  motions,  finite  depth,  irregular  fluid  bottom  geometry, 
axisymmetric  geometry,  and  fluid  stratification.   Bai  also 
provided  criteria  for  properly  placing  the  radiation  boundary 
and  developed  a  new  form  of  the  radiation  condition  for  the 
axi-symmetric  case. 

The  author  has  extended  the  studies  of  Chenault 
and  Bai  in  the  solution  of  the  first-order  problems  with  an 
improved  version  of  Chenault 's  computer  program.   The  results 
are  presented  later  in  this  work  along  with  results  separately 
reported  by  Newton,  Chenault,  and  Smith  [23]. 
2.   Second  Order  Solutions 

a.   General  Potential  Theory 

The  origins  of  the  second-order  theory  of  gravity 
waves  begin  again  with  the  work  of  Stokes  and  his  classic 
solution  of  nonlinear  traveling  waves.   Others  have  extended 
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the  work  of  Stokes,  but  it  appears  that  a  modern  version 
of  Stokes'  theory  was  first  compiled  by  Stoker  [35]  in  1957. 
The  generally  accepted  terminology  associated  with  Stoker's 
work  is  "perturbation  theory"  which  of  course  has  applications 
in  many  fields . 

Fontannet  [5],  in  1961,  appears  to  be  the  first 
to  attempt  to  use  the  basic  nonlinear  perturbation  theory  as 
described  by  Stoker  to  solve  the  problem  of  finding  the 
characteristics  of  waves  generated  by  a  plane  wave-maker  using 
Lagrangian  coordinates.   Ogilvie  [25],  in  1963,  found  the 
first  and  second-order  forces  on  a  fixed  or  oscillating 
circular  cylinder  submerged  in  a  fluid  with  a  free  surface. 
However,  his  second-order  forces  were  only  those  due  to  the 
first-order  solution,  and  only  the  time-independent  part  was 
reported.   Tuck  [37]  found  the  second-order  forces  on  a 
submerged  cylinder  in  a  uniform  stream  in  1965.   Salvesen  [33] 
extended  Tuck's  work  to  include  wii^g-shaped  bodies  in  his 
Ph.D.  dissertation  in  1966. 

Lee  [17]  and  Parissis  [27]  appear  to  be  the  first 
to  apply  second-order  perturbation  theory  to  attempt  to  solve 
the  fluid-structure  interaction  potential  flow  problem  with 
a  free  surface  in  1967.   Lee  confined  his  studies  to  the 
heaving  of  floating  circular  and  U-shaped  sections,  while 
Parissis  studied  only  the  circular  shape  in  heave;  both 
studies  were  in  fluids  of  infinite  extent.   These  investiga- 
tions, as  well  as  those  previously  mentioned,  required  the 
development  of  a  fluid-structure  second-order  boundary 
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condition  (discussed  in  more  detail  later)  in  addition  to 
the  free-surface  equations  compiled  by  Stoker. 

The  work  of  Lee  and  Parissis  was  followed  by  that 
of  Potash  [31]  in  1970.   Potash  extended  the  second-order 
studies  to  include  the  additional  degrees  of  freedom  of  roll 
and  sway  and  the  attendant  coupling  between  modes  of  motion. 
Recently,  Garrison  [7]  completed  a  second-order  solution  of 
the  problem  of  determining  the  dynamic  loads  experienced  by 
a  "fixed"  three-dimensional  body  due  to  an  incoming  wave 
train.   His  theory  includes  submerged  and  surface-piercing 
bodies . 

All  of  the  second-order  solutions  previously 
mentioned  used  some  form  of  singularity  distribution  in  the 
flow  field.   The  functional  forms  of  the  singular  functions 
used  were  usually  taken  from  the  work  of  Wehausen  and  Laitone 
[41]. 

b.   Finite  Element  Solutions 

The  author  has  been  able  to  find  only  one 
solution  to  the  second-order,  of  the  class  of  problems  being 
discussed,  by  the  finite  element  method.   Allouard  and  Coudert 
[2]  presented  a  paper  at  the  International  Symposium  on  Finite 
Element  Methods  in  Flow  Problems  in  January  of  1974.   They 
presented  a  method  of  solving  problems  related  to  floating 
bodies  in  water  waves  under  unsteady-linear  or  stationary- 
nonlinear  (second-order)  conditions. 
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B.   OBJECTIVES  OF  PRESENT  RESEARCH 

The  experimental  studies  of  Paulling  and  Porter  [29], 
Paulling  and  Richardson  [28],  and  Vugts  [40],  previously 
mentioned,  indicated  good  agreement  with  the  developed  linear 
theory  if  the  amplitude  of  motion  was  small.   Potash  [31] 
observed  that  small  variations  in  amplitude  caused  large 
discrepancies  between  Vugt ' s  experimental  results  and  those 
of  linear  theory.   Potash  [31]  has  shown  that  second-order 
effects  become  significant  at  higher  frequencies  and  are  more 
significant  when  roll  and/or  sway  modes  are  present.   However, 
he  did  not  examine  the  effect  of  finite  depth  on  second-order 
solutions . 

The  present  work  has  two  primary  goals.   The  first  is  to 
determine  if  the  finite  element  method  can  be  successfully 
applied  to  finding  solutions  to  second-order  boundary  value 
problems  of  the  type  previously  described.   And,  if  success 
is  achieved  in  the  first  objective,  the  second  major  goal  is 
to  study  the  significance  of  nonlinear  effects  in  heave  in 
water  of  finite  depth. 
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II.   FORMULATION  OF  THE  MATHEMATICAL  MODEL 

The  physical  problem  is  to  determine  the  force  required 
to  sustain  steady,  vertical  harmonic  motion  of  a  floating 
cylindric  ship  hull  form.   The  ship  will  be  considered 
oscillating  in  an  "ideal"  fluid  of  infinite  horizontal  extent 
and  arbitrary  constant  depth. 
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FIGURE  1.   Hull  at  Free  Surface 

A.   THE  NONLINEAR  BOUNDARY  VALUE  PROBLEM 

The  domain  of  the  boundary  value  problem  is  shown  in 
figure  1.   In  the  figure,  b  is  the  ship's  half  beam,  h  is  the 
mean  fluid  depth,  d  is  the  ship's  draft,  and  2A  is  the  total 
submerged  cross-sectional  area  of  the  ship. 

The  fluid  is  assumed  to  be  incompressible,  inviscid, 
homogeneous,  and  without  surface  tension.   It  is  well  known 
that  under  the  assumption  of  irrotationality ,  a  velocity 
potential  exists.   Furthermore,  the  potential  will  be  a 
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single-valued  function  in  any  simply-connected  region.   The 
convention  for  the  potential  used  herein  is 

$$  =  ut  +  vj  ,  (1) 

where  u  is  the  x-component  of  fluid  velocity,  v  is  the 
y-component  of  fluid  velocity,  and  1  and  j  are  the  unit 
vectors  in  the  positive  right-handed  cartesian  coordinate 
system  directions  x  and  y  respectively.   The  symbol  V  is  the 
standard  "gradient"  used  in  vector  calculus . 

The  nature  of  the  problem  dictates  that  at  a  large 
horizontal  distance  from  the  ship's  hull  (cross-hatched  area 
fig.  1)  only  an  outgoing  one-dimensional  traveling  wave  will 
be  present.   Therefore,  it  is  advantageous  to  consider  the 
reduced  region  shown  in  figure  2.   In  figure  2,  w  is  the 
region  semi-width.   The  meanings  of  the  other  symbols  shown 
in  figure  2  will  become  apparent  in  the  following  discussion 
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FIGURE  2.   Reduced  Region 
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Conservation  of  mass  requires: 

V2$  =  0  ,    in  R   .  (2) 

p 

R  is  the  fluid  domain  and  V  is  the  Laplacian  operator. 

There  are  two  nonlinear  free  surface  conditions,  one  dynamic 
and  one  kinematic.   The  dynamic  condition,  from  the  Bernoulli 
equation,  is 


•t  +  \     [*£  +     S>p     +  E  +  yg    =    0   .  (3)1 


Equation  3  must  hold  on  the  moving  free  surface  (p  =  0) .   In 
the  equation,  g  is  the  acceleration  of  gravity,  p  is  the 
overpressure  (referred  to  atmospheric),  p  is  the  fluid 
density,  and  t  is  real  time.   The  subscripts  denote  partial 
differentiation  with  respect  to  the  indicated  variable.   The 
kinematic  condition  comes  from  a  basic  assumption  from 
continuum  mechanics  simply  stated  by  Stoker:   "A  particle 
once  on  the  free  surface  remains  on  it."   This  assumption 
translates  into  the  following  expression 


D[n(x,t)  -  y]   Q  m 

Dt         U   *  v  ; 


In  equation  3,  the  possible  pure  function  of  time, 
which  may  appear,  is  assumed  to  be  identically  zero  without 
loss  of  generality  (See  Stoker  [35],  p.  10). 
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Again,  equation  *!  must  be  satisfied  on  the  moving  free 
surface.   In  equation  4,  D/Dt  denotes  the  co-moving  derivative 
(sometimes  called  the  substantial  or  total  derivative)  and  the 
function  n(x,t)  is  defined  by  the  following  expression 

y  =  n(x,t)  .  (5) 

Therefore,  n  is  defined  to  be  the  wave  height  above  the  line 
y  =  0  at  any  point  x  and  any  time  t. 

Equation  4  can  be  rewritten  using  the  definition  of  $  in 

p 
equation  1  and  the  definition  of  the  co-moving  derivative 

on  the  moving  free  surface . 

On  and  "near"  the  bounding  surface  Sk  (radiation  boundary) 

<Kx,y,t)  -*  $*(x,y,t)  ,   as  x+oo   .        (7) 

In  equation  7,  $*  denotes  the  velocity  potential  of  a 
traveling  wave  moving  away  from  the  ship.   It  must  be  mentioned 
here  that  inherent  in  the  statement  of  equation  7  is  the 
assumption  that  the  surface  Sk  is  "far  enough"  away  from  the 
ship's  hull.   This  point  will  be  discussed  in  more  detail 
later. 

Turning  to  the  interface  Sp  between  the  fluid  and  the 
ship's  hull,  the  kinematic  condition  is 


2 

The  co-moving  derivative  Is  defined  as 

D(  )/Dt  =  u(  )x  +  v(  )   +  (  )t   . 
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;$-n  =  v-n   .  (8) 


The  vector  n  is  the  unit  outward  normal  vector  (from  the 
region  R)  to  the  surface  at  the  point  of  interest .   The  dot 
is  the  standard  scalar  product  and  V  represents  the  vector 
velocity  of  the  ship  at  any  instant  of  time. 

Since  the  motion  of  the  ship  is  assumed  to  be  restricted 
to  the  vertical  direction  (heave),  the  bounding  surface  S- 
becomes  a  plane  of  symmetry.   The  surface  S,-,  represents  the 
bottom  of  the  fluid  region  and  is  rigid.   Therefore,  both 
surfaces  have  the  same  kinematic  boundary  condition 

$$.n  =  0  .  (9) 

An  examination  of  equations  3  and  6  indicates  the  nonlineari- 
ties  inherent  in  the  boundary  value  problem.   An  additional 
problem,  not  as  obvious,  is  the  fact  that  the  free  surface 
and  the  hull  are  moving  boundaries.   In  its  present  form  the 
problem  has  not  yielded  solutions.   Therefore,  some  approxi- 
mate theory  is  applied. 

In  the  development  which  follows,  a  perturbation  analysis 
will  be  applied  which  when  carried  to  the  second-order  will 
produce  two  linear  boundary  value  problems  from  the  one 
presently  posed.   Further,  the  boundary  conditions  of  both 
linear  problems  will  be  referred  to  fixed  boundaries .   In 
addition,  appropriate  equations  for  determining  force,  and 
wave  profiles  will  be  obtained. 
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B.   PERTURBATION  ANALYSIS 

First  it  is  assumed  that  $  and  n  each  possess  a  perturba- 
tion series  representation  as  follows 

4>  =  e$(1)  +  eV2)  +  0(g3)  ,  (10) 

n  =  en(1)  +  eV2)  +  0(e3)  ,  (11) 

where  z   is  a  small  non-dimensional  parameter  which  is  a 
measure  of  the  size  of  the  ship's  motion  and  is  defined  as 
e  =  a/b  ("a"  is  the  amplitude  of  the  ship's  motion  and  "b" 
is  the  ship's  half-beam).   Also  in  equation  10,  $    and  $ 
are  defined  to  be  the  first  and  second-order  velocity 
potentials,  respectively,  and  similarly  for  r\  and  r\         . 

Note  that  as  e-*0  there  is  no  disturbance  and  therefore  there 
are  no  zero-order  functions.   Substitution  of  equation  10 
into  equation  2  yields 

eV2$(1)  +  e2vV2)  +  0(e3)  =  0  .  (12) 

Equation  12  implies  that 


V2$(1)  =  0  ,   and  (13) 

V2<D(2)  =  0  ,   in  R  .  (14) 


A  similar  result  holds  for  equation  9. 
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1 .   The  Linearized  Free  Surface  Condition 

The  following  development  of  the  linearized  free 
surface  conditions  follows  in  general  that  of  Stoker  [35]. 
The  final  linearized  equations  for  the  free  surface  boundary 
condition  are  obtained  in  the  following  steps.   Equations  5, 
10,  and  11  are  substituted  into  equation  3,  noting  that  p  =  0 
on  the  free  surface.   Also,  it  is  observed  that  $/:   ,  $    , 
and  $    may  be  represented  by  the  following  form  of  the 
Taylor  series 


$v(x,y,t)  =  $v(x,0,t)  +  n(x,t)$   (x,0,t) 

2 

+  \   <Dvyy(x,o,t)  +  o(n3)   .  (15)3 


After  collecting  coefficients  of  like  powers  of  e  in  the 
resulting  equation,  the  following  set  of  equations  results 


gn(1)   +  <^1}     =     o     ,  (16) 


gn(2)    +    #(2)    +n(D#(l)    +|  [{$(D}2+    {f(l)}2]    =    0    f 

on  S3    .       (17) 

Applying  a  very  similar  procedure  to  equation  H   yields  the 
following  result  to  the  second  order 


The  subscript  v  denotes  partial  differentiation  with 
respect  to  x,  y,  or  t. 
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n(t1)  -  ^X)      -  0   ,  and,  (18) 


n<2)  _  #(2)  +  n(l),(l)  .  n(D,(l)   =  o  ,   on  5,  .   (19) 
't      y      x   x         yy        "      3 


Differentiating  equation  16  with  respect  to  time,  and 
eliminating  n.  from  equations  16  and  18  yields  the  well 
known  linearized  free  surface  condition 

g$y1)  +  *tt*  =  °  '  on  y  =  °  •         (20) 

In  equation  20,  because  of  the  use  of  the  Taylor  series 
expansions  in  n,  the  functions  are  evaluated  for  y  =  0, 
rather  than  on  the  actual  free  surface. 
Following  the  same  procedure,  as  before,  for  equations 
17  and  19  yields 

g*(2)  +  ,(2)  =  1  #(D(#(1)  +  g$(D) 

S  y      tt    g  vt   ^tty   s  yy  ; 

-  2(^(1)^,1)  +  t^U^h       .  (21) 

x   xt     y   yt 

Note  the  nonlinear  function  of  $    on  the  right-hand  side  of 
equation  21. 

2.   The  Ship  Interface  Condition 

Since  the  ship  boundary  is  also  in  "motion",  a 
similar  procedure  must  be  applied  there  to  "fix"  the  boundary 
Sp.   The  development  follows  to  some  extent  that  of  Potash 
[31].   Consider  figure  3  in  which  the  general  ship  hull  form 
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FIGURE  3.   Hull  Coordinate  System 


is  described  in  two-dimensional  cartesian  coordinates  in 
terms  of  the  parameter  s  (arc  length)  in  the  following  way 


x(s) 


y  =  y(s)  . 


(22) 


The  inertial  coordinate  system  is  the  same  as  shown  in 
figures  1  and  2.   The  parametrically  represented  body 
coordinates  (equation  22)  correspond  to  the  inertial  coordi- 
nates of  the  body  when  the  ship  is  at  rest  (neutrally  buoyant 
position) .   The  motion  of  the  ship  in  this  reference  frame 
can  be  defined  in  the  following  way 


x(s,t)  =  x(s)  , 


y(s,t)  =  y(s)  +  e  yh(t) 


(23) 
(24) 
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where  the  function  y,  is  defined  to  be 

h 

yh  =  Re{belat}  .  (25) 

In  equation  25,  a  is  the  circular  frequency  of  the  harmonic 
motion  and  1  is  the  imaginary  unit. 

The  velocity  vector  for  the  ship  referred  to  the 
inertial  reference  frame  is 

V  =  e  yh  j  .  (26) 

The  dot  denotes  ordinary  time  differentiation.   The  unit 
outward  normal  vector  (from  the  region)  in  terms  of  the 
parametric  coordinates  becomes 

n  =  -y'l  +  x'j  ,  (27) 

where  the  prime  indicates  differentiation  with  respect  to 
the  variable  s.   Substituting  equations  26  and  27  into 
equation  7  yields  the  ship-fluid  interface  boundary  condition 

?*«n  =  eyhx'  .  (28) 

Now  expanding  $  in  a  Taylor  series  in  e  (similar  to  free 
surface  expansion)  yields 

*[x(s,t),y(s,t),t]  =  $(x,y,t)  +  eyh$y(x,y,t)  +  0(e2)  .   (29) 

Using  equations  10  and  29,  the  normal  derivative  of  $  on  S2 
may  be  expressed  as 
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$$«n  =  (-y'<l>(1'  +  x'<i>(1))e 
x        y 

+  [-yf^2)  +  ^,$y2)  -  y'y^xy5  +  7'yh$yy)]6:2  +  0(e3)   » 

on  S2  .      (30) 

Substituting  equation  30  into  equation  28  and  equating  like 
powers  of  e  yields  the  following  set  of  equations  valid  on 
the  fixed  boundary  S?: 

V<D(1)  -n  =  yhx"'   ,  (3D 

v*(2U  =  -x'yh«<J>  +  y'yh^\  (32) 

Equation  32  may  also  be  written 


V*(2)-n  =  -yh(h.V$(1))y   .  (33) 


It  remains  to  consider  the  boundary  S^  (radiation 
boundary)  and  the  function  $*.   For  consistency,  $*  must 
also  be  expanded  in  a  perturbation  series 

<i>*  =  e$*(1)  +  e2$*(2)  +  0(e3)  .  (3*0 

Potash  [31]  and  others  have  shown  that  $*    and  $*    ,  for 
a  region  of  infinite  depth,  each  represent  a  simple  harmonic 
traveling  wave  of  appropriate  frequency  and  wave  length. 
Zienkiewicz  and  Newton  [43]  have  provided  the  appropriate 
homogeneous  boundary  condition  for  the  surface  Sh  to  success- 
fully pass  a  wave  form  without  reflection.   The  equations  are 
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$M<D.S  =  -  J-  ♦.(!>    ,  (35) 


Cl 


^«(2).S  =   -  i  ♦•(25    .  (36) 


C2 


In  the  above  equations  c-,  and  c?  (the  wave  celerities)  are 
defined  by  the  following  implicit  relationships 


c,  =  (g/a)  tanh  (2&)   ,  (37) 

■la  {s   -. 


cp  =  (g/2a)tanh(^-)   .  (38) 


3.   Perturbation  Formulas  for  Pressure, 
Force  and  Wave  Amplitude 

The  determination  of  dynamic  loads  on  the  ship 
requires  the  following  equations  which  follow  from  the  pertur- 
bation series  assumed  for  $  and  n 

p(x,y,t)  =  p(0)(x,y)  +  ep(1)(x,y,t)  +  e2p( 2) (x,y ,t )  +  0(e3)  , 

(39) 
F(t)  =  eF(1)(t)  +  e2F(2)(t)  +  0(e3)  .  (40) 

In  equation  40,  F(t)  is  defined  as  the  total  force  per  unit 
length  required  to  sustain  simple  harmonic  vertical  trans- 
lation (see  figure  3). 

Following  a  similar  procedure  as  before,  evaluating 
the  Bernoulli  equation  (equation  3)  on  the  ship's  hull  using 
equations  23,  24,  and  29  yields 
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p[x(s),y(s),t]  =  gy(s)  +  e[gyh(t)  +  •f1*(x,y,t)] 


1 
P 


+  e2[yh^}  +  *£2)  +  |  ({^1}}2  +  {^1}}2)]  +  0(e3)  .  (41) 


By  comparing  equation  39  with  equation  41  we  find 


P(0)  =  -  Pgy(s)  ,  (42) 


P(1)  =  -  P[gyh  +  *£1}]  ,  (43) 

P(2)  --  P[yh*<J>  +  |({^}2  +  C^}2)  +  *<2>]  . 

(44) 

At  this  point  observe  that  the  second-order  pressure  is  made 
up  of  contributions  from  both  the  first  and  second-order 
potential  solutions. 

Potash  [31]  and  others  have  provided  relations  for 
the  force  per  unit  length  acting  on  the  cylinder.   However, 
in  this  treatment,  it  is  desired  to  obtain  relations  for  the 
exciting  force  as  previously  defined.   The  motive  is  that  the 
theoretical  result  obtained  should  correspond  to  that 
measured  by  an  appropriate  physical  experiment.   The  relation 
between  the  exciting  force  F,  harmonic  displacement  y,  and 
hydrodynamic  pressure  p,  is 


F(t)  =  ep2Ay,  -   /  pn-j  ds  ,  (45) 

-I 
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where  I   is  one-half  of  the  total  symmetric  arc-length  of  the 
ship's  hull  (see  figure  3).   It  should  be  noted  here  that 
Potash  [31]  has  examined  the  question  of  the  time-dependent 
total  arc-length  of  the  actual  wetted  hull  and  its  effect  on 
the  limits  of  integration  in  equation  45.   However,  It  will 
be  stated  without  proof  that,  provided  the  sides  of  the  ship 
are  vertical  at  y  =  0,  the  second-order  effect  is  zero.   The 
only  hull  form  examined  to  the  second  order  satisfies  this 
stated  requirement.   Therefore,  equation  45  is  exact  for  the 
case  considered. 

In  the  calculation  of  the  force  F,  the  constant 
hydrostatic  lifting  force  which  comes  from  p    will  not  be 
considered  since  it  is  not  in  the  definition  of  F.   Accordingly, 
the  following  expression  for  F(t)  results  using  equations  27, 
39,  43,  44,  and  45 

F(t)  =  e[p2Ayh  -   /   -p(gyh  +  ^1))x«  ds  ] 

~ ■  At 

+  e2[-P   /  -(yn^y}  +  l[Ux1)}2  +  {*y1)}^  +  $t2))I'  ds]* 
~*  (46) 

Comparing  equation  46  with  equation  40  yields  the  following 
set  of  equations 

F(1)(t)  =  P2Ayh  +  p  /  (gyh  +  ^1})x'  ds   ,  (47) 

F(2)(t)  -  P  Avg}  + 1«41}>2  +  <*<1}>2)  +  ^2)]x«  ds  . 

(48) 
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The  wave  amplitude  may  be  written  by  inspection  of 
equations  16  and  17. 

n(1)(x,t)  =  -  I  ^1}  ,  (H9) 

n<2>(*,t)  =  -  I[*<2>  +  n(1)^J)  +  |({»<1))2  +  (^1}>2)]  . 

(50) 

The  boundary  value  problems  for  $    and  $    may  now 
be  formally  stated.   However,  at  this  point,  a  standard 
separation  of  variables  technique  is  applied  to  eliminate 
time  from  the  solution  (since  the  solution  is  assumed 
periodic) .   The  appropriate  definitions  are 


*(1)(x,y,t)  =  ReU(1)(x,y)  elat}  ,  (51) 


*(2)(x,y,t)  =  ReU(2)(x,y)  eiat}  .  (52) 


Note  the  newly  defined  spatial  potential  functions  <}>    and 

(2) 

<f>    are  in  general  complex.   Also  because  of  the  definitions 

in  equations  51  and  52  it  is  necessary  to  replace  y,  by  the 


complex  form 


y*  =  beiat   .  (53) 
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C.  THE  BOUNDARY  VALUE  PROBLEM  FOR  <t>(1) 

The  boundary  value  for  <fr  '  is  formally  stated  using 
equations  9,  13,  20,  25,  31,  and  35. 

V2<f>(1)  =  0  ,   in  R   ,  (5l|) 

g4»y1}    -    aV1}    =   0    ,      on  S3  (55) 

V"<J>(1).n   =   iabx'     ,      on  S2  (56) 

^(l).S  =  i2|^l     ,     onS,  (57) 

Cl  H 

V<f>(1).n   =   0    ,      onS1   and  S5    .  (58) 

D.  THE  BOUNDARY  VALUE  PROBLEM  FOR  4>^2) 

(2) 

Similarly,  the  boundary  value  problem  for  4>  may  be 

formally  stated  using  equations  9,  1^,  21,  25,  32,  and  36. 

V2(j>(2)  =  0  ,   in  R  ,  (59) 

-  ia[{0>^1)}2  +  {^1)}2]  ,   on  S2  (60)" 

V2>-n  =  ^-^  ,   o„S2  (61) 


The  factor  H   on  the  right-hand  side  of  equations  60  and 
6l  is  due  to  complex  algebra  involved  and  is  explained  in 
Appendix  A. 


39 


V2>-S  =  -2i|iHi  ,   onS,  (62) 

?V2  -n  =  0  ,   onS1  and  S5  (63) 


Note  the  striking  similarities  between  the  first  and  second- 
order  problems.   Note  also,  that  the  nonlinearities  in  the 
original  problem  have  been  approximated  by  the  nonlinear 
terms  in  <J>    in  equation  60,  and  the  nonlinear  terms  in  y, 

and  cj)    in  equation  6l  (y,  does  not  explicitly  appear)  . 

2 
The  additional  nonlinear  amplitude  dependence  comes  from  e 

in  the  original  perturbation  expansions  (equations  10  and 
11) .   The  structure  of  the  problem  requires  the  solution  of 
<J>    first  in  order  to  define  the  solution  for  <J> 

The  conversion  to  complex  algebra  allows  another  observa- 
tion to  be  made.   Since  only  the  real  part  of  <j>    is 
meaningful,  Appendix  A  demonstrates  the  required  complex 

manipulation.   Close  examination  of  that  result  indicates 

(2) 
that  the  nonhomogeneous  -boundary  conditions  on  <p    contain 

(2) 

time-independent  terms.   This  implies  <t>    has  a  time- 
dependent  and  a  time-independent  part.   However,  since  the 

physical  quantities  of  interest  in  the  solution  all  involve 

(2)  (2) 

derivatives  of  f   ,  the  time-independent  solution  of  <J> 

5 
will  not  be  pursued  further.    The  above  comments  are  also 

true  for  the  second-order  functions  of  pressure,  force,  and 


Lee  [18]  comments  on  the  time-independent  solution  in 
terms  of  mass  transport. 


HO 


wave  amplitude.   However,  the  time-independent  quantities 

are  of  interest  and  have  physical  meaning  as  will  be  seen 

later. 

The  boundary  value  problem  for  <pK        is  valid  only  in 
i 
water  of  infinite  depth.   The  following  discussion  will 

explain  the  problem  and  develop  a  boundary  value  problem 

(2) 
which  will  allow  for  the  solution  of  <Jr    in  finite  depth. 

E.   THE  SECOND-ORDER  BOUNDARY  VALUE  PROBLEM 
WITH  FINITE  DEPTH 

Potash  [31]  has  shown  that  the  asymptotic  form  of  the 

velocity  potential  $  as  x -*  °°  becomes  (in  water  of  infinite 

depth) 

f(-.fl.t)  =  Re{EG1eiC°t-(klx^(1)"+  ^^^t-i^^y^hl  > 

(64) 

where  k..  denotes  the  wave  number  of  the  first  order  wave 
(k1  =  a/c-,  )  and  G-,  and  G2  are  coefficients  which  depend  on 
frequency  a.   However  in  the  more  general  case  of  finite 
depth,  an  additional  term  due  to  the  Stokes'   second-order 
potential  appears  as  follows 

*/   a  *.\        u    f    u   i[ot-(kix+Y    )]  L  2ru   -i(kpX+Y    ) 
<K°°,.0,t)  =  Re{eH-.e       -1   '    ,J   +  e  [Hpe    £ 

-  3  _3£^  H12cosh(2k1h)     2i(k1x+Y(1))n   2iat,   .    (65) 
8c^     sinh2(k1h)   e  J  e    ' 


The  Stokes  potential  to  the  second  order  is  well  known 
to  be  zero  in  infinite  depth. 
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In  equation  65,  k~  is  defined  as  the  second-order  wave  number 
(kp  =  2a/c2).   H,  and  Hp  are  constants  similar  to  G-,  and  G? . 
Close  examination  of  equation  65  reveals  that  the  radiation 
boundary  condition  (equation  62)  would  fail  since  the 
"Stokes  wave"  which  is  present  is  traveling  at  the  same 
celerity  as  the  first-order  solution.   This  means  that  the 
radiation  boundary  would  "reflect"  the  Stokes*  wave  since 
the  boundary  condition  defined  by  equation  62  will  only  pass 
a  wave  traveling  at  a  celerity  of  Cp,  defined  by  equation  38. 
A  successful  resolution  of  the  above  described  difficulty 

requires  the  formulation  of  a  different  boundary  value 

(2) 

problem.   We  define  the  complex  velocity  potential  i>  as 

t(2)  .  .  «£_  H12cosh[2k1(y+h)]  ^(j^U))  ^      (^&) 

Qc±2  sinh2(k1h) 


which  corresponds  to  the  third  term  in  equation  65.   Then, 


following  previous  conventions,  the  complex  velocity  potential 

_(2) 

<|)    is  defined  as 

*<*>  -  ♦<*>  -  *<2>   .  (66) 


1.   The  Second-Order  Boundary  Value 
__ 

Problem  for  <j> ;  J 

(2) 

The  definition  of  "$"    in  equation  66  allows  the 


formal  statement  of  the  following  boundary  value  problem, 

the 
(2) 


(2) 
which  follows  from  equation  66,  the  definition  of  ^    and 


the  boundary  value  problem  for  <}> 
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V2   ,f>(2)    =0      ,      in  R    ,  (67) 


bj<2>  -  toV2)  =  i  ^  [■«\!1)  *  8il!»] 


2g  Yy  &yyy 


-  iali^h2   +  {♦J1'}2]  -  g*(2)  ♦  4aV2)  ,  on  83, 


(68) 


btlV     ,„        b*<!> 


V2).S-f  (^S-+^))  _x'(^+^2))  ,  onS2> 

(69) 

*<2>.J  -  _  2i£i^  ,  ons4  ,  (70) 

+_(2)  ■>     ±  (2)  •> 

Vcj)^;-n  =  -  Vij/^-n   ,   onS.   ,  (71) 


^cf)(2)  -n  =  0  ,   on  S^   .  (72) 


Equation  70  will  now  properly  behave  because  the  asymptotic 

—(2) 
form  of  (j)v    using  equations  65  and  66  is 

T<2>(.,y)  e2i0t  -  H2  e-2i(k2x+Y(2))  e2iat  (?3) 


Note  the  celerity  of  the  complex  wave  form  ~§        e     is  c~ 
and  the  radiation  boundary  condition  will  now  work. 


F.   FORMULAS  FOR  PRESSURE,  FORCE,  AND  WAVE  AMPLITUDE 

The  development  which  follows  requires  additional  defini- 
tions, necessary  because  of  the  separation  of  variables  scheme 
previously  introduced.   They  are: 
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p(1)(x(y,t)  =  Re{f(1)elot}  ,  p(2>  =  Re{p<2)e2iot} ,  (7*) 
F(1)(t)  =  Re{f(1)elat}  ,  F(2)(t)  =  Re{f (2)e2iat }  ,  (75) 
n(1)(x,t)  =  Re<n(1>(x)elot>  ,   n(2)  =  Re{n<2)(x)e21ot}.  (76) 


The  functions  p^   ,  p    ,  rf   ,  n"    are  complex.   This  is 
also  true  for  the  constants  f    and  f 

Using  the  above  equations  and  equations  43,  44,  47,  48, 
49,  and  50  the  equations  for  pressure,  force,  and  wave 
height  for  both  first  and  second-order  solutions  follow. 

p(1)  =  -p[bg  +  ia<{>(1)]  ,  .      (77) 

p(2)  .  _p[^4^  +  ^(♦<1)>2  +  {^1}}2)  +  2icV2>]  ,    (78) 

f(1)  =  -2pbAa2  +  p   /  (gb  +  lo^1')**    ds  ,  (79) 

-I 


f(?)  =  _p   /(l0*y1)b  +  kuilh2  +  i^h2i  +  2io^2')x-  dS> 

(80) 


~2 '  2Llvx  •  l*y 


-(1)  =  _  lof(1)  (81) 

g 

1       f2l    ian(1)d)(1)    1  {*x1)}2    {*v1)}2 
-(2)  =  -  i  (2ia<f>(2)  +  iqn    fy    +  k— ^ +  — 5L ])  . 

(82) 
Again,  equations  78,  80,  and  82  include  the  special  considera- 
tion of  the  complex  algebra  involved  which  is  explained  in 
Appendix  A. 
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III.   FINITE  ELEMENT  DISCRETIZATION 

The  Finite  Element  Method  and  its  application  to  scalar 
field  problems  is  described  in  various  texts,  including  one 
by  Zienkiewicz  [45].   The  essential  ideas  involved  follow. 
The  region  of  figure  2  is  replaced  by  a  grid  or  mesh  of 
finite  elements  as  shown  in  figure  k. 


FIGURE  4.   Finite  Element  Mesh 

Then  it  is  assumed  that  the  field  function  sought  can  be 
represented  by  the  following  equation 

4>(x,y)  =  N  J   ,   in  R  .  (83) 

The  total  number  of  nodes  in  the  region  will  be  defined  to 
be  m.   Then  N  is  a  lxm  row  vector  of  shape  (interpolating) 
functions  and  <J>  is  an  m  x  1  column  vector  of  nodal  values  of 
the  complex  potential  $  (in  this  problem) .   Equation  83 
implies  that  knowledge  of  the  vector  <J>  defines  <j>(x,y) 
everywhere  in  R. 
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There  are  several  choices  of  shape  functions  N.  which  may 
be  used  in  equation  83.   This  work  selected  the  cubic  shape 
functions  described  by  Zienkiewicz  [*J5]  as  the  "serendipity 
family".   The  elements  used  are  isoparametric.   A  brief 
description  of  the  12-noded  isoparametric  elements  is  given 
in  Appendix  B  for  the  unfamiliar  reader.   The  key  advantage 
of  the  isoparametric  element  lies  in  the  fact  that  regions 
with  curved  boundaries  may  be  readily  represented. 

The  development  that  follows  defines  the  set  of  linear 
equations  which  determine  the  vector  <J>  of  equation  83. 

One  general  approach  to  discretization  of  several  types 
of  field  problems  in  finite  element  analysis  is  to  define  a 
functional  whose  integral  over  the  domain  of  the  problem  is 
to  be  minimized.   The  Calculus  of  Variations  shows  that  if  the 
Euler  equation  of  the  integral  to  be  minimized  is  the  governing 

equation  (i.e.,  Laplace's  Equation)  then  the  two  problems  are 

7 
equivalent. 

However  Zienkiewicz  and  Newton  [43]  have  shown  the 

discretization  may  be  readily  accomplished  by  applying  the 

Galerkin  process.   The  weight  functions  are  chosen  to  be  the 

nodal  shape  functions.   Mathematically  stated 


/  NTUXX  +  d>   )  dR  =   0  .  (84) 

R 


This  statement  assumes  the  same  boundary  conditions  are 
used. 
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The  superscript  T  denotes  transposition  and  the  function  <J> 
is  considered  to  be  any  of  the  complex  potential  functions 
whose  solution  is  sought.   Application  of  the  divergence 
theorem  transforms  equation  8*1  into 


'  [!?x  fyl 
R    A   y 


dR   =   /  N  V<J)-n  dS   , 
S  ~ 


(85) 


where  S  denotes  the  boundary  of  the  region  R. 

A.   THE  DISCRETIZED  FIRST-ORDER  PROBLEM 

Application  of  equations  84  and  85  to  the  first-order 
boundary  value  problem  (<j)  )  stated  in  Chapter  II  yields 
the  following  matrix  equation 


(_02Q  +  i£  D  +  H)  *(1)  =  iabr(1)  . 


(86) 


1  = 


The  definitions  of  the  matrices  in  equation  86  follow.   The 
matrix  H  is  m x m,  real,  symmetric,  and  defined  as 


H  = 


/  [NT  NT] 
R   ~X  ~y 


N 
~x 

N„ 


dR 


(87) 


The  matrix  H  comes  from  the  integral  on  the  left-hand  side 
of  equation  85.   The  matrix  D  is  mxm,  real,  symmetric,  and 
defined  by 


D  =   /    NTN  dS  . 
z  S„   ~  ~ 


(88) 
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D  is  contributed  by  the  boundary  integral  on  the  right-hand 
side  of  equation  85  and  represents  the  discretization  of  the 
radiation  boundary  condition  (equation  57) . 

The  matrix  QQ  is  real,  -m x m,  symmetric  and  given  by 


QQ  =  I  /   NTN  dS  .  (89) 

z  s3 

Q  parallels  D  in  its  origin  and  represents  the  contribution 

~0  as 

from  the  homogeneous  free  surface  boundary  condition  (equation 
55). 


The  vector  r    is  real  and  defined  by 


(1)  =  r        mT-, 
2 


vKX)    =   f        N  x»  dS  .  (90) 

S. 


The  above  vector  is  the  discretized  ship-fluid  interface 
kinematic  boundary  condition  (equation  56).   The  surfaces 
S,  ,  S~,  S|,  and  S,-  have  homogeneous  boundary  conditions  and 
therefore  make  no  contribution  to  r    .   It  follows  that  r 
is  non-zero  only  at  nodes  along  the  ship's  hull  (surface  Sp)  . 

B.   THE  DISCRETIZED  SECOND-ORDER  PROBLEM 

-(2) 
The  problem  for  <j> v    only  will  be  presented  since  it 

( 2) 

contains  the  solution  for  <pK    '    in  the  limit.   The  boundary 

—(2) 
value  problem  for  <fr   yields  a  linear  system  of  complex 

equations  similar  to  equation  85. 


(-Ho2Q„   +  2ia  D  +  H)  -(2)  =  r(2)   ^  (91) 

^O      Cq~     2~         ~ 


H8  ( 


(2) 
The  right-hand  side  vector  r    in  equation  80  receives 

contributions  from  the  surfaces  S,  ,  S~  and  S^.      Accordingly 

a  further  decomposition  is  necessary: 


^(2)  .  ^(2)  +  rU)    +  ^(2)  _  (92) 

(2) 

The  vectors  r^   are  defined  below.   First  to  be  consistent, 

(2)  —(2) 

the    function   ty  must   be   represented  in  the   same  way   as   <$> 


*(2)   =  N  i>{2)  (93) 


( 2) 
Then  v^.        becomes,  from  equations  71  and  85, 


r{2)    =  /  NTN  ij,(2)  dS  .  (94) 

^1 

(2) 

The  nodal  values  of  ip    are  determined  by  partial  differen- 

(2) 
tiation  of  the  function  ty  and  then  representing  the 

(2) 

function  \p  in  similar  manner  as  equation  93. 

(2 ) 
The  vector  v\        is  defined  by  equations  69  and  85 

cj>(1)b  <f>(1)b 

r<2)  =  /   NTNy"'[^—  +  ij/2)]dS  -  /   NTNx'  l=M —  +  ^2)]dS  . 

2  2  (95) 

Similarly  r~    '    follows  from  equations  68  and  85 

j(2)  =  I  /   nTN[-  ^  a  +  *  g  -  lay  -  g^2)  +  4c2^2)]  dS  . 
3  (96) 
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The  i   component  of  a  is 


i   Yi   yi,y 


(97) 


where  <J>.    and  <f>.    are  the  i   components  of  <Jr  ^  and  <J>    , 
l        l ,y  ~        ~y 

respectively.   Similarly 


Hi       yi       yi,yy 


and 


(98) 


Tl  ■  <>>a  ♦  <>>2 


(99) 


The  nodal  values  of  the  various  partial  derivatives  of  <f> 
are  determined  by  differentiating  equation  83.   The  method 
of  calculation  is  discussed  in  the  next  section.   The  reader, 
unfamiliar  with  isoparametric  elements,  should  refer  to 
Appendix  B  prior  to  reading  the  next  section. 
1.   Calculation  of  Nodal  Values  of  the 


■  ■ 


Field  Derivatives  of  ft 


OT 


The  vectors  ft^1  \    K      >    *xv  >  and  ^vv   reQuired  in 
equations  95  and  96   are  determined  by  applying  standard 
concepts  of  finite  element  analysis.   The  appropriate  relation 
for  the  i   node  of  a  given  element  is 


* 


4> 


(l)e 

i,x 

(De 


=  J 


-1 


■rf 

Ne 

(De 


(100) 
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J  is  the  2x2,  real  Jacobian  matrix  of  the  transformation. 
The  vector  N  is  a  1  x 12  row  vector  of  element  level  shape 
functions.   The  superscript  "e"  denotes  all  quantities  in 

brackets  are  at  the  element  level.   In  equation  100,  the 

-1      e 
elements  of  J   and  N   are  evaluated  at  node  i. 

Application  of  equation  100  to  all  twelve  nodes  of  a 

(l)e      (l)e 
given  element  determines  the  vectors  $  and  <}> 

~x       ~y 

(De     (De 
Then,  replacing  <$>  by  <J>^    in  equation  100  yields  the 

(De      (l)e 
components  of  the  vectors  <f>     and  <f>Jtv    .   Mathematically 

stated 


<$> 


4> 


(De 

i,xy 

(De 

i,yy 


r-i 


Ne 


(De 


(101) 
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IV.   THE  COMPUTER  PROGRAM  AND  NUMERICAL  SOLUTION 

A.   MAIN  FEATURES  OF  THE  COMPUTER  PROGRAM 

The  formation  and  solution  of  the  linear  system  of  equations 
defined  by  equation  86  ( <jr    system)  or  equation  91  (<J> 
system)  was  accomplished  on  an  IBM  360/67  computer  at  the  Naval 
Postgraduate  School.   The  computer  has  a  maximum  word  length 
of  sixteen  bytes. 

A  comparison  of  equations  86  and  91  demonstrates  the  close 
similarity  of  the  coefficient  matrices  for  the  first  and 
second-order  solutions.   Note  that  only  the  coefficients  of 
Q  and  D  change.   An  additional  observation  follows  from 
examination  of  equations  87,  88,  and  89  defining  H,  D,  and 
Q  respectively.   Note  that  all  three  matrices  are  purely  a 
function  of  mesh  geometry.   Therefore,  once  a  mesh  is  chosen 
to  define  the  region  of  interest,  the  three  matrices  previously 
mentioned  need  only  be  calculated  once  regardless  of  the  forced 
motion . 

A  well  known  property  of  the  finite  element  method  is  that 
it  usually  produces  a  banded,  symmetric,  coefficient  matrix. 
That  is  the  case  here  and  this  has  obvious  computational 
advantages.   However,  in  this  problem,  the  final  system  of 
equations  is  a  function  of  the  frequency  of  motion  as  for 
example  in  the  first-order  problem 


K^V1'  =  iobr(1)   .  (102) 
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The  matrix  K    is  an  m x m  complex  matrix  defined  by  equation 
86  to  be 


K(1)  =  -a2Q^  +  —  D  +  H   .  (103) 

z  z°        ci  ~    z 

Although  K^    is  complex  and  frequency  dependent,  it  is 
possible  to  perform  most  of  the  triangular  decomposition  of 
K   '  (by  Gauss  elimination)  only  once,  using  real  arithmetic, 
for  a  given  problem  geometry.   This  economization  was  first 
employed  by  Chenault  [4]  and  is  described  in  more  detail  in 
Appendix  C.   In  addition,  this  author  employed  a  vector 
storage  of  the  upper  triangular  portion  of  K    ,  which  greatly 
reduced  the  storage  requirements  of  the  program.   The  vector 
storage  scheme  made  it  possible  to  store  only  the  portions  of 
K   '  which  are  non-zero. 

~  o 

The  net  advantage  of  employing  Chenault ' s   scheme  and  the 
compact  vector  storage  was  that  the  system  of  733  complex- 
equations  which  is  generated  by  the  mesh  shown  in  figure  4 
could  be  solved  for  both  the  first  and  second-order  problem 
twenty  times  in  twenty  minutes.   Further,  only  ten  percent  of 
the  total  matrix  was  actually  stored  and  the  computer  solution 
was  accomplished  in  core. 


o 

The  Chenault  scheme  described  in  Appendix  C  was  developed 
by  Professor  R.  E.  Newton  at  the  Naval  Postgraduate  School. 
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B.   NUMERICAL  SOLUTIONS 

The  solution  of  the  first-order  problem  for  a  given  mesh 
geometry  is  required  before  the  second-order  solution  can  be 
obtained.   Further,  the  first-order  solution  must  be  highly 
accurate  if  a  good  second-order  solution  is  to  be  expected. 
This  fact  is  made  more  obvious  by  re-examining  equations  68 
and  69  in  part  II . 

1.   First  Order  Numerical  Solutions 

Several  first-order  problems  were  solved  and  closely 
examined  for  accuracy.   The  numerical  solutions  are  presented 
in  part  V. 

a.   Mesh  Requirements 

One  of  the  primary  considerations  in  choosing  a 
proper  mesh  concerns  the  positioning  of  the  radiation 
boundary  (Sj,).   It  may  be  recalled  that  the  surface  Sh  must 
be  far  enough  away  from  the  ship  hull  (location  of  disturbance) 
so  that  the  local  effect  of  the  disturbance  has  decayed  and 
only  an  outgoing  wave  is  present.   Bai  [3]  has  addressed  this 
question  analytically.   He  has  shown  that  the  eigenfunctions 
associated  with  the  non-propagating  velocity  field  decay 
exponentially  with  x.   Bai  recommends  the  criterion 

w  "    =  ah  +  b  ,  (104) 

min  * 

where  a  is  a  monotonically  increasing  function  of  frequency. 
The  range  of  a  is  from  1.5  (low  frequency)  to  3.   The  values 
of  a  were  determined  based  on  the  decay  rate  of  the  most 
persistent  eigenfunction  and  an  attenuation  factor  of  0.01. 

5^ 


This  investigation  has  found  the  above  criterion  somewhat 
conservative.   A  value  of  a  =  2/3  was  found  to  be  adequate 
for  the  velocity  potential  to  numerically  approach  its 
asymptotic  limit  for  all  frequencies  studied  for  the  semi- 
circular hull  (h/b  =  6).   Similarly,  a  value  of  a  =  9/10  was 
found  to  be  adequate  for  the  bulb  hull  form  described  in 
part  V  (h/b  =  10)  for  all  frequencies  studied.   The  values 
of  a  specified  above  were  determined  for  the  deepest  regions 
studied  for  each  hull  form.   The  region  widths  used  in  this 
case  (w/b  =  6,  w/b  =10)  were  found  to  be  more  than  adequate 
for  shallower  depths . 

Another  requirement  on  w  which  tends  to  restrict 
its  maximum  possible  value  comes  from  the  practical  limits  on 
mesh  fineness  and  the  resulting  computer  storage  requirements. 
The  mesh  must  be  sufficiently  fine  to  represent  properly  a 
traveling  surface  wave.   The  finer  the  mesh,  the  larger  the 
computer  storage  requirements  and  the  longer  the  computer 
solution  time.   Therefore,  this  latter  restriction  tends  to 
oppose  the  former  one. 

Bai  [31  used  quadratic  elements  and  recommended 
at  least  10  surface  nodes  per  wave  length.   Visser  and  Van 
der  Wilt  [39]  also  used  quadratic  elements  and  recommended 
that  an  element  could  not  span  more  than  one  quarter  of  a 
wave  length.   This  investigation  examined  this  question  in 
detail.   As  previously  mentioned,  cubic  isoparametric  elements 
were  used.   It  was  found  that  if  the  gravity  wave  did  not  have 
to  be  propagated  more  than  two  wave  lengths,  the  elements 
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could  be  allowed  to  span  six-tenths  of  a  wave  length. 
However,  as  in  the  case  Visser  and  Van  der  Wilt  [39]  studied, 
the  successful  propagation  of  a  gravity  wave  over  many  wave 
lengths  requires  a  finer  mesh  spacing.   The  one  quarter  wave 
length  Visser  and  Van  der  Wilt  recommend  was  found  to  be  a 
practical  upper  limit. 

A  parallel  problem  to  positioning  the  radiation 
boundary  concerns  the  location  of  the  bottom  (surface  S._) 
to  simulate  infinite  depth.   Comparison  of  the  finite  element 
solutions  obtained  in  this  study  with  analytic  solutions 
based  on  infinite  depth  led  to  the  conclusion  that  a  keel  to 
bottom  clearance  (h  -  d)  of  5b  will  adequately  simulate 
infinite  depth.   However,  satisfying  the  classical  hydro- 
dynamic  requirement  that  the  depth  must  be  at  least  one-half 
of  a  wave  length  becomes  impractical  at  low  frequencies. 
Ignoring  this  requirement  did  not  seem  to  affect  solution 
accuracy. 

b.   Mesh  Data  Preparation 

The  meshes  for  a  given  boundary  geometry  were 
prepared  using  a  very  versatile  mesh  generator  computer 
program  developed  by  Adamek  [1]  at  the  Naval  Postgraduate 
School.   This  tool  greatly  expedited  mesh  preparation.   For 
example,  the  entire  data  for  the  mesh  shown  in  figure  H 
(733  nodes,  132  elements)  could  be  generated  in  twenty 
seconds.   A  modification  was  used  to  supply  accurate  values 
of  the  nodal  coordinates  along  the  ship-fluid  interface. 
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c.   Convergence  of  First-Order  Solutions 

Convergence  of  the  first-order  solutions  could 
not  be  determined  by  standard  numerical  techniques  (i.e., 
mesh  refinement)  due  to  computer  size  limitations.   However, 
the  excellent  agreement  between  the  first-order  solutions 
obtained  in  this  study  and  both  theory  and  experiment  left  no 
doubt  that  the  solutions  were  quite  good. 
2.   Second-Order  Numerical  Solutions 

The  successful  solution  of  the  first-order  problem 
made  study  of  the  second-order  problem  possible. 

a.   Mesh  Requirements 

All  of  the  requirements  on  mesh  geometry  previously 
mentioned  for  the  first-order  solutions  of  course  apply  to  the 
second-order  solution.   However,  the  second-order  problem 
demands  even  more  stringent  requirements.   One:  reason  for  this 
follows  from  the  fact  that  the  wave  lengths  of  the  propagating 
portion  of  the  second-order  velocity  potential  are  .much 
shorter  (one-fourth  of  the  first-order  wave  lengths  in  infinite 
depth  solutions).   Therefore,  since  it  is  impractical  to  use 
a  different  mesh  for  each  solution,  the  mesh  must  have  a  free 
surface  element  spacing  approximately  one-fourth  of  that 
required  for  the  solution  of  <j>    for  a  given  frequency  range 
to  be  studied.   The  mesh  shown  in  figure  k   was  used  to  solve 
the  infinite  depth  second-order  problem.   It  contains  132 
elements  and  733  nodes. 
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b.   The  Second-Order  Boundary  Condition  on  S_ 

Additional  constraints  on  mesh  geometry  were 
found  to  be  required.   These  constraints  arise  from  the 
requirements  to  calculate  second-order  partial  derivatives  of 
<J)    which  appear  in  equations  95  and  98.   It  was  essential 
to  have  rectangular  elements  with  equal  node  spacing  along 
the  free  surface  for  success  in  representing  the  non- 
homogeneous  boundary  condition  given  by  equation  96. 

The  further  discussion  of  this  point  will  be 
aided  by  the  definitions  of  some  complex  functions.   Define 

Q(x,o)  =  isgl  r-c2^  +  g$>]  -  icru^)2  +  u^)2:  . 

(105) 
R(x,0)    =    g^2)    -    4oV2)       .  -  (106) 

In  terms  of  these  definitions  the  right-hand  side  of  equation 
68  is 

P(x,0)  =  Q(x,0)  -  R(x,0)   .  (107) 

The  reader  should  recall  that  the  successful  solution  of  the 

(2) 

boundary  value  problem  for  <J>    using  the  method  of  finite 

elements  required  the  definition  of  a  new  boundary  value 

—(2) 
problem  for  the  function  <$>K       .   The  reason  was  that  in  the 


finite  depth  problem  a  second-order  Stokes'  wave  is  produced 

,S  W' 

■(2) 


with  celerity  c,  as  well  as  a  second-order  wave  with  celerity 


Cp  (contained  in  IJr   ).   The  presence  of  the  two  waves  with 
different  celerities  meant  the  radiation  boundary  condition 
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would  fail.   The  problem  was  mathematically  eliminated  by 

subtracting  the  well-known  analytic  solution  for  the  Stokes 

(2) 
wave  portion  of  the  second-order  solution  of  <p      '    on  every 

boundary  of  the  region  R.   The  mathematical  steps  are  straight- 
forward as  shown  in  Chapter  II.   The  asymptotic  value  of  the 
function  P(x,0)  is  zero  as  x  gets  large.   But,  the  function 
R(x,0)  which  represents  the  classic  Stokes  second-order  wave 
free  surface  boundary  condition  has  the  form 

R(X)0)  =    310*    |t(l)(Wj0)|2  e-2i(k1X+Y(1))  >  (108) 
g  sinh  k-,d 

everywhere  on  S~.  Therefore,  if  the  numerical  representation 
of  the  function  P(x,0)  is  to  approach  zero  near  the  radiation 
boundary;  the  numerical  representation  of  the  function  Q(x,0) 
must  approach  R(x,0)  or 

Q(x,0)  -*-  R(x,0)    as   x  ■*■  °°  (109) 

Note  that  the  amplitude  and  phase  of  <f>    must  be  known  to 

(2) 

completely  define  R(x,0)  and  therefore  ijr   . 

The  above  discussion  points  out  the  necessity  of 

a  very  accurate  numerical  determination  of  the  function  Q(x,0). 

—  (2) 
The  second-order  solutions  for  <f>    were  closely  examined  near 

the  radiation  boundary  and  it  was  found  that  on  the  average 

the  function  Q(x,0)  was  within  one-percent  of  the  analytic 

function  R(x,0).   This  accuracy  was  achieved  even  though 

formation  of  Q(x,0)  requires  calculation  of  second  derivatives 
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of  <f>    using  element  shape  functions  which  guarantee  only 

C°  continuity  across  inter-element  boundaries. 

c.   The  Ship-Fluid  Interface  Second-Order 
Boundary  Condition 

The  accurate  representation  of  the  ship-fluid 
interface  (S?)  boundary  condition  for  the  second-order 
boundary  value  problem  (equation  69)  is  of  course  essential. 
Potash  [31]  encountered  difficulties  in  making  a  numerical 
calculation  on  this  boundary  because  the  potential  singulari- 
ties which  define  his  solution  lie  along  S? .   This  results 
in  the  required  functions  4>    and  <J>    being  only  piecewise 
continuous.   This  problem  required  an  additional  approximation 
to  be  made. 

The  method  used  herein  also  has  the  problem  that 

<f>    ,  <P    >  <p    j  and  <J>    are  not  continuous  across  element 
Tx   '  Yy   »  Yxy  '     Tyy 

boundaries.   However,  no  singularities  are  present.   Extensive 
investigations  in  this  work  have  shown  that  it  is  essential 
that  the  elements  along  the  hull  be  as  square  as  possible, 
and  further  that  the  node  spacing  on  the  elements  along  the 
hull  should  be  uniform.   It  was  found  that  if  the  hull 
elements  were  distorted  (not  square  or  rectangular)  the 
numerical  representation  of  the  functions  <J>    and  <b    was 
highly  unstable.   However,  the  mesh  requirements  for 
numerically  calculating  <j>    and  <K    along  the  hull  were 
much  less  stringent. 

The  adequacy  of  the  hull  elements  used  in  this 
investigation  (see  figure  4)  was  tested  by  defining  a  known 
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potential  function  on  the  elements  along  S? 

♦*  =  |  y2  +  i  xy  (110) 

and  then  applying  the  computer  program's  numerical  scheme 

t       t 
to  calculate  the  functions  <J>   and  $   .   Observe  that 


*L  =  i   >  (in) 


and 


xy 


4=1.  (112) 


It  was  found  that  the  cubic,  isoparametric  elements  along  S? 

produced  an  average  error  of  about  one  percent  in  calculating 

the  functions  <J>   and  <\>      .   The  maximum  error  observed  was 

xy     Tyy 

five  percent  at  isolated  nodes. 

d.   Convergence  of  Second-Order  Solutions 

Again,  computer  capacity  limitations  prevented 
the  application  of  systematic  mesh  refinement  to  determine 
convergence  of  the  second-order  solutions.   However,  it  was 
possible  to  refine  the  elements  next  to  the  ship's  hull  and 
no  significant  change  in  the  solutions  was  observed  (the 
elements  along  the  free  surface  were  performing  properly) . 
Further,  it  is  shown  in  Chapter  V  that  second-order  results 
for  infinite  depth  are  in  satisfactory  agreement  with  those 
of  other  researchers.   No  comparison  is  available  for  finite 
depth  solutions,  but  a  smooth  transition  is  observed  as  the 
depth  becomes  effectively  infinite  (See  Chapter  V) . 
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V.   NUMERICAL  RESULTS 

A.   DEFINITIONS 

The  definitions  which  follow  are  useful  for  presentation 
of  results . 

1.   Added  Mass  and  Damping  Coefficients 

The  complex  amplitude  F^    of  the  first-order 
resultant  hydrodynamic  force  acting  vertically  on  the  ship's 
hull  is 

pd)  =  /  p^x'  ds  .  (113) 

y     -I 

The  real  part  of  F    is  in  phase  with  ship's  acceleration 
and  the  ship's  added  mass  is  defined  by 

Re{F^1}} 
ma  =    /     ,  (114) 

c  a 

2 
where  the  ship's  acceleration  amplitude  is  given  by  a  a. 

The  added  mass  coefficient  is  defined  to  be  the  ratio  of 

added  mass  to  ship's  displaced  mass 

m     Re{F(1)} 
m   2pA    2pAaa2 

The  imaginary  part  of  F^   is  in  phase  with  the 

Q 

velocity  of  the  ship  and  is  the  "non-conservative"^ 


The  potential  theory  is,  of  course,  conservative.  However, 
the  waves  which  propagate  to  "infinity"  represent  work  done  by 
the  ship  which  is  effectively  lost. 
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component.   Accordingly,  the  damping  coefficient  is  defined 
to  be 

-ImCF^1*) 

Cd  =  1 .  (116) 

2pAacT 

C  and  C,  are  frequency  dependent  for  a  given  hull  form  and 
region  depth. 

A  dimensionless  frequency  parameter  is  defined  as 

follows 

2. 
6  -  2g2.  .  (117) 

The  parameter  6  is  commonly  used  by  other  researchers  and  is 

used  herein  for  ease  of  comparison  of  solutions .   The 

parameter  6  is  equivalent  to  —j—   in  infinite  depth  (A  is  the 

wave  length) . 

The  added  mass  and  damping  coefficients  previously 

defined  are  for  heave.   The  same  coefficients  for  sway  motion 

are  similarly  defined.   The  only  difference  in  the  definitions 

is  F    is  replaced  by  the  horizontal  force  amplitude  F 
y  x 

defined  by 

F(1)  =  -  /   p(1)y~'  ds  ,  (118) 

x       -I 

where  p^  '    is  the  hydrodynamic  pressure  due  to  the  swaying 
motion.   There  is  a  coupling  between  sway  and  rolling  motion 
for  a  general  ship  cross-section.   However,  the  only  hull 
form  studied  in  the  sway  mode  was  semi-circular  in 
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cross-section  and  the  coupling  force  for  this  case  is  zero. 
Therefore,  the  appropriate  coefficients  typically  used  are 
not  defined  herein. 

One  further  comment  is  warranted  here.   The  second- 
order  time-dependent  hydrodynamic  forces  are  harmonic  with 
frequency  2a.   Therefore  the  definition  of  second-order  added 
mass  and  damping  coefficients  is  meaningless. 
2.   Dimensionless  Force  Coefficients 

The  force  coefficients  f    and  f    are  also 

presented  in  dimensionless  form.   The  non-dimensionalizing 

p 
factor  for  the  various  force  amplitudes  is  1/pgb  .   Table  I 

presents  the  correspondence  between  the  dimensional  and  non- 
dimensional  force  quantities.   The  second-order  quantities 

(2)       (2) 

f    and  fj,   which  appear  in  table  I  are  amplitudes  of  the 
s        d 

time-independent  (static)  and  time-dependent  (dynamic) 

(2)  (2) 

portions  of  f    respectively.   The  existence  of  f^    is 

readily  verified  by  examining  equations  75  and  80  in  light  of 

the  discussion  in  Appendix  A. 


TABLE  I 
DIMENSIONLESS  FORCE  COEFFICIENTS 


Force  Amplitude        f(1)       f£2)       f^2) 


Dimensionless  F(1)       F^2)       F.(d: 


Force  Amplitude 


d 
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3 .   Wave  Amplitude  at  Infinity 

Another  suitable  measure  of  the  damping  of  ship's 
motion  (to  the  first-order)  is  the  wave  amplitude  at 
"infinity"  (i.e.,  near  the  radiation  boundary).   Further, 
experimentally  it  has  been  found  to  be  an  easier  quantity 
to  measure  (particularly  for  high  frequency  motion  since  the 
damping  force  is  very  small  then).   Accordingly,  the  dimen- 
sionless  first-order  wave  amplitude  at  infinity  is  defined 
by 

-(l>_]^_!™  (ll9) 

B.   HULL  FORMS 

Other  investigators  have  studied  a  variety  of  hull  forms 
both  by  experiment  and  analytically.   The  semi-circular  hull 
was'  chosen  because  of  the  wealth  of  data  available  for 
"comparison  of  both  first  and  second-order  solutions.   The 
bulb-shaped  hull  studied  by  Paulling  and  Porter  [29]  was 
also  examined  for  first-order  solutions.   The  shape  of  the 
bulb  hull  in  comparison  to  the  semi-circular  hull  is  shown 
in  figure  5. 

The  coordinates  of  the  bulb  hull  studied  are  defined  by 
the  following  mapping 


e    e    e 

x  +  iy"=^  +  -^+-4+-|   •  (120) 
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Centerline 


Waterline 


Semi-Circle,  d/b  - 1 


FIGURE  5.   Hull  Forms 

The  complex  number  r,   represents  the  points  on  a  circle  in 
the  £-plane.   The  choice  of  radius  in  the  c; -plane  defines 
the  scale  of  the  hull  cross-section.   The  values  of  the  e. 
are  given  in  table  II. 


TABLE 

II 

CONSTANTS  FOR 

HULL  FORMS 

Hull 

el 

e3 

e5 

Semi-Circular 

0.0 

0.0 

0.0 

Bulb 

-0.7132 

-0.02096 

0.0 

A/bd 

0.785^ 
0.6950 
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C.   FIRST-ORDER  NUMERICAL  SOLUTIONS 
1.   Semi-Circular  Hull 

The  semi-circular  hull  form  has  been  studied  by  many 

researchers  mentioned  previously  in  part  I.   Specifically, 

the  original  solution  by  Ursell  [38]  for  heaving  motion  has 

been  confirmed  analytically  by  Porter  [30],  Vugts  [40],  and 

others.   Ursell's  solution  has  also  been  experimentally 

confirmed  by  Paulling  and  Porter  [29]  and  Vugts.   The 

analytic  results  of  Ursell  for  C  and  C-,  are  shown  by  the 
J  m      d  J 

continuous  curve  in  figure  6.   The  numerical  results  obtained 
by  the  finite  element  method  (PEM)  in  this  work  are  repre- 
sented by  the  circled  points  in  figure  6.   The  FEM  points 
shown  were  obtained  using  a  mesh  similar  to  that  shown  in 
figure  4,  but  using  only  62  cubic  elements  and  359  nodes. 

The  extension  of  the  solutions  for  C   and  C,  to 

m      d 

finite  depth  is  easily  accomplished  by  simply  redefining  the 
mesh  geometry  (decreasing  the  vertical  dimension) .   Kim  [14] 
and  Yu  and  Ursell  [42]  reported  analytic  results.   The 
results  obtained  by  the  FEM  in  this  work  for  h/d  =  4  agree 
with  those  reported  by  Kim  for  C   and  n*^   within  three 
percent  over  the  range  0.25  <  6  <  5.   The  results  of  Kim  are 
in  disagreement  with  those  obtained  by  Yu  and  Ursell  for 
h/d  =  2.   The  results  obtained  herein  fall  between  the  two. 
Bai  [3]  reports  inability  to  verify  Kim's  result,  but  gives 
no  details.   Bai  used  the  finite  element  method  to  obtain 
his  solutions  as  mentioned  in  part  I. 
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FIGURE  6.   Added  Mass  and  Damping  Coefficients  in 

Heave,  Semi-Circular  Hull,  Infinite  Depth 


Figure  7  shows  the  effect  of  finite  depth  on  C 
obtained  in  this  work.   The  figure  shows  that  significant 
changes  do  not  occur  until  very  shallow  water  is  encountered 
(h/d  <_   2)  .   Further,  the  effect  of  shallow  water  on  C   is 
more  pronounced  at  very  low  and  very  high  values  of  6 .   Note 
the  unusual  variation  with  depth  for  6=0.5.   C  decreases 
with  decreasing  depth  (opposite  the  trend  at  higher  values 


68 


\ 


1.25 


1.00 


0.75 


m 


0.50 


0.25 


0 


s 

4.5- 

3.0  iL  \^ 

2.510^.    ^^4^_ 

20rC^^r? 1 

l-5T^^^_^ 

l.oi _, 

J 

j 

0.5<i 

4 
h/d 


8 


FIGURE  7.   Variation  of  Added  Mass  Coefficient  with 
Depth,  Semi-Circular  Hull  in  Heave 
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of  6).   Although  it  does  not  appear  true  in  figure  7,  C 
(for  <5  =  0.5)  has  reached  its  asymptotic  value  at  h/d  =  6. 
Figure  8  shows  the  corresponding  effect  of  finite  depth  on 
C,.   The  figure  indicates  that  C,  is  only  sensitive  to 
shallow  water  effects  at  very  low  values  of  6 . 
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FIGURE  8.   Variation  of  Damping  Coefficient  with 
Depth,  Semi-Circular  Hull  in  Heave 
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The  semi-circular  hull  was  also  studied  in  sway 
motion  to  the  first-order.   The  only  changes  in  the  boundary 
value  problem  for  <jr    defined  in  part  II  are  in  equation  56 
which  defines  the  surface  S, .   The  changes  are 

V<J)(1).n  =  iabx'  ,   on  S2   ,  (121) 

for  equation  56,  and 

<J>(1)  =  0   ,   on  S1   ,  (122) 

for  equation  58.   The  boundary  condition  on  S,-  defined  by 
equation  58  remains  unchanged.   The  ship's  motion  is  defined 
in  a  manner  analogous  to  equations  23  and  24 


x(s,t)  =  x(s)  +  e  xh(t)  ,  (123) 

y(s,t)  =  y(s)  .  (124) 

The  function  x,  (t)  is  defined  the  same  as  yh(t) 

iat 
xh  =  Reibe1  z]    ,  (125) 

and  e  is  defined  in  the  case  of  sway  motion  to  be 


a 
e  =  b 


(126) 


The  symbol  "a"  is  the  amplitude  of  the  sway  motion  in  this 

case . 

Vugts  [40]  presents  both  experimental  and  analytic 

solutions  for  C  and  C,  for  the  semi-circular  hull  in  sway. 

m      d 

Figure  9  shows  a  comparison  of  Vugts'  theoretical  solution 
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FIGURE  9.   Added  Mass  and  Damping  Coefficients  in  Sway, 
Semi-Circular  Hull,  Infinite  Depth 
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with  the  FEM  results  obtained  herein.   The  FEM  calculated 
points  are  shown  by  the  solid  circles  and  squares  in  the 
figure.   The  solid  curve  (from  Vugts'  theory)  was  scaled 
from  a  very  small  graph  and  may  be  somewhat  in  error. 

2.  Bulb  Hull 

Paulling  and  Porter  [29]  reported  both  analytic  and 
experimental  results  for  the  bulb  hull  form.   They  show 
excellent  agreement  by  comparing  first-order  exciting  force 
amplitude  coefficients  F^    and  wave  amplitude  coefficients 
rj"^   respectively.   Figures  10  and  11  compare  the  values  of 
F^  '    and  rp    obtained  by  the  FEM  herein  (circled  points) 
and  Paulling  and  Porter's  theoretical  calculations  (solid 
line).   The  FEM  mesh  parameters  were  h/d  =  2  and  w/b  =  25. 
This  value  of  w/b  was  found  later  to  be  quite  conservative 
(w/b  =  10  would  have  been  adequate) . 

Figure  12  shows  added  mass  and  damping  coefficients 
for  the  bulb  hull  vs  h/d  for  various  values  of  6.   The  curves 

have  a  similar  trend  to  those  of  the  semi-circular  hull  form. 

-(1) 
The  corresponding  curves  of  force  Fv    and  wave  amplitude 

n^   are  not  presented.   However,  the  maximum  change  in 

exciting  force  amplitude  does  not  exceed  five  percent  for 

h/d  =  1.2  as  compared  to  infinite  depth  (h/d  =  2).   Similarly, 

the  maximum  change  in  rf^   was  only  ten  percent  over  the  same 

range  of  h/d. 

3.  Wave  Generator  Studies 

MacCamy  [20],  in  1957,  reported  an  analytic  theory 
using  a  Green's  function  approach  for  solving  the  problems 
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FIGURE  10.   Exciting  Force  Amplitude  in  Heave, 
Bulb  Hull,  Infinite  Depth 
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FIGURE  11.   Wave  Amplitude  for  Heave, 
Bulb  Hull,  Infinite  Depth 
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FIGURE  12.   Variation  of  Added  Mass  and  Damping 

Coefficients  with  Depth,  Bulb  Hull  in  Heave 


of  gravity  wave  production  by  a  moving  partition  and  gravity 
wave  reflection  from  a  horizontal  strip.   The  problem  of 
gravity  wave  production  by  a  moving  partition  was  studied  in 
this  work  to  further  help  determine  proper  mesh  characteris- 
tics for  gravity  wave  propagation. 

The  wave  generation  problem  is  readily  placed  in  the 
framework  of  the  class  of  problems  being  studied  herein  with 
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only  minor  modifications.  °  All  of  the  boundary  conditions 
posed  for  <J>    in  part  II  remain  unchanged  except  those  on 
S1   and  S2.   The  surfaces  S-[  and  S2  become  a  vertical  line. 
The  region  therefore  is  rectangular.   The  left-hand  boundary 
will  be  designated  as  S'  as  shown  in  figure  13. 
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FIGURE  13.   Reduced  Region  for  Wave  Generator  Problem 


The  boundary  condition  on  Si  becomes 


v>(1)-n  =  q(y)  ,   on  Sj_   , 


(127) 


The  function  q(y)  describes  the  mode  of  displacement  of  the 
wave  generator.   The  plunger  type  and  flapper  type  (shown  in 
figure  13)  wave  generators  were  studied  and  compared  with 


This  observation  is  not  true  in  the  potential  theory 
singularity  distribution  approach.   The  source  singularity 
must  be  replaced  by  doublets  and  the  problem  completely 
reformulated. 
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MacCamy's  results.   The  function  q  (y)  for  the  plunger  type 
generator  is  simply 


qp(y)  =  U0  ,  (128) 

where  y   is  the  maximum  amplitude  of  the  wave  generator 
motion.   Similarly,  the  function  qf(y)  for  the  flapper  type 
generator  is 


qf(y)  =  ^j-21  yQ  .  (129) 


The  comparison  of  MacCamy's  results  with  those  obtained 
herein  is  shown  in  figure  1H    for  the  plunger.   The  circled 
points  were  obtained  by  the  FEM  and  the  solid  line  is  MacCamy's 
theoretical  curve.   Similarly,  the  comparison  of  MacCamy's 
result  for  the  flapper  type  wave  generator  vs  the  FEM  results 
obtained  in  this  work  is  shown  in  figure  15.   The  excellent 
^agreement  allows  two  observations  to  be  made.   First,  the  FEM 
approach  to  solution  of  problems  of  the  type  being  presented 
is  highly  versatile  in  its  application.   Second,  the  radiation 
boundary  placement  was  further  verified. 

D.   SECOND-ORDER  NUMERICAL  RESULTS 

1.   Semi-Circular  Hull  —  Infinite  Depth 

The  first-order  solutions  just  presented  established 
confidence  that  second-order  solutions  could  be  attempted. 
The  only  hull  form  studied,  to  the  second-order,  was  semi- 
circular.  The  accuracy  of  the  second-order  solutions  was 
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FIGURE  lH.      Wave  Amplitude  Generated  by  a 
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FIGURE  15.   Wave  Amplitude  Generated  by  a 
Flapper  Type  Wave  Generator 
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established  by  comparing  with  the  work  of  other  researchers. 
Lee  [17,18],  Parissis  [27],  and  Potash  [31]  have  reported 

numerical  solutions  for  second-order  exciting  force  coeffi- 

_(2)   —(2)  (2) 

cients  Fd   ,  Fg    and  phase  angle  y   .  Figure  16  compares  the 

—(2)  —(2) 

values  of  F^   and  F^  '    reported  by  Lee,  Parissis,  and 

Potash  with  the  numerical  results  obtained  by  the  FEM  in  this 

work.   The  points  presented  from  the  other  researchers'  works 

were  all  obtained  from  Potash  [31].   Potash  indicates  Lee's 

points  are  corrected  in  some  manner  from  his  original  work 

[17]  and  otherwise  unpublished. 

Figure  16  shows  excellent  agreement  for  F^    and  very 

5 

—(2) 
good  agreement  for  F^   .   Only  Potash's  points  show  signi- 
ficant deviation  and  then  only  for  6  near  1.5.   Potash 
reported  some  numerical  difficulties  for  a  range  of  6  near  1 
and  therefore  Potash's  results  are  not  presented  for  those 

values  of  6  which  are  suspected  to  be  in  error.   Similarly, 

( 2 ) 
figure  17  shows  a  comparison  of  y  (phase  angle)  results. 

The  figure  shows  the  agreement  is  not  as  good.   However,  none 

of  the  other  works  shown  (Lee,  Parissis,  Potash)  agree  in  any 

consistent  fashion  over  a  large  range  of  6.   The  reason  for 

the  discrepancy  has  not  been  resolved. 

2.   Semi-Circular  Hull  —  Finite  Depth 

—(2) 
Figure  18  presents  the  effect  of  depth  on  F^   and 

~(  2) 

Fv  '    as  obtained  by  the  FEM  in  this  work.   The  author  knows 
s 

of  no  other  solutions.   Solutions  were  obtained  for  values 
of  h/d  from  6.0  (infinite  depth)  to  1.5.   The  figure  shows 
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FIGURE  16.   Second-Order  Exciting  Force  Amplitude  in 
Heave,  Semi-Circular  Hull,  Infinite  Depth 
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FIGURE  17.   Second-Order  Exciting  Force  Phase  Angle 

in  Heave,  Semi-Circular  Hull,  Infinite  Depth 
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FIGURE  18.   Variation  of  Second-Order  Exciting 
Force  Amplitudes  with  Depth, 
Semi-Circular  Hull  in  Heave 
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— (2)  —(2) 

that  Fl   and  Fv  J    are  most  sensitive  to  finite  depth  at  high 
d        s 


values  of  6  or  in  very  shallow  water  (h/d  <  1.5). 


(2) 


Figure  19  shows  the  corresponding  effect  on  y  of 

finite  depth  for  the  same  range  of  h/d.   Note  that  at 

(2) 
h/d  =  l.M  6  has  negligible  effect  on  y         . 

Also  the  effect  of  decreasing  depth  causes  a  monotonic 

(2) 

variation  of  y  at  fixed  6. 


FIGURE  19.   Variation  of  Second-Order  Exciting  Force  Phase 
Angles  with  Depth,  Semi-Circular  Hull  in  Heave 
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Figure  20  presents  the  ratio  of  F^2'  to  f'1' 
plotted  vs  6  for  various  values  of  h/d.   The  curves  peak 
close  to  6  =  1  for  all  values  of  h/d  studied.   The  curves 
peak  because  the  first-order  force  F^  '  passes  through  a 
minimum  very  near  6=1  for  all  values  of  h/d  studied. 


FIGURE  20. 


Ratio  of  Second-Order  Exciting  Force  Amplitude 
to  First-Order  Exciting  Force  Amplitude, 
Semi-Circular  Hull  in  Heave 
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A  quantitative  illustration  of  the  physical  signif- 
icance of  the  nonlinear  effects  (second-order)  is  readily 
made  using  figure  20.   Assume  a  ship  with  an  approximately 
circular  hull  form  is  being  excited  such  that  the  maximum 
acceleration  of  the  motion  is  0.2  g's.   The  maximum  time- 
dependent  amplitude  of  the  force  is  given  by 


F(t)|     =  e|f(1)|  +  e2|f(2)|   ,  (130) 

i max     i     i      i     i   > 


which  may  be  written 


F(t)l  py  w(2) 

'max  m   1  +  £  P.  ^  (131) 


-77m--  iTb  JflT  ' 

Assume  h/d  =  1.5  and  6  =  1.0,  the  maximum  ship  acceleration 
(in  g's)  may  be  written 

HiLtl)]  =  e  6  .  (132) 

g 

Equation  132  implies  e  =  0.2  under  the  assumed  conditions. 
Using  figure  20  and  equations  131  and  132 


P(t)  I 

jj2££   =  1  +  (0.2X1.75)  =  1.35  (133) 


e  f 


Therefore,  neglecting  the  second-order  contribution  could 
result  in  a  thirty-five  percent  error  in  lF(t)lmax  • 

Figure  21  further  demonstrates  the  nonlinear  (second- 
order)  effects  on  F(t).   The  function  F(t)/ef(1)  is  plotted 
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vs  t.   The  ship's  displacement  as  a  function  of  time  is  also 
shown  on  the  plot  for  reference  (cosine  function  of  unit 
amplitude).   The  parameters  chosen  for  the  plot  are  h/d  =  1.5, 
6  =  1.0  and  e  =  0.2  (corresponds  to  an  equivalent  acceleration 
of  0.2  g's) . 
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FIGURE  21.   Total  Exciting  Force  vs.  Time 
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Figure  22  illustrates,  for  the  same  values  of  the 
parameters,  the  nonlinear  effects  on  the  gravity  wave 
generated  by  the  ship's  disturbance.   The  figure  presents 
the  sum  of  the  dimensionless  first-order  wave  amplitude  and 
the  second-order  Stokes'  wave  (both  have  celerity  c, ) . 
Separately,  the  second-order  wave  from  "<p  '    (celerity  cp)  is 
plotted.   Both  curves  are  plotted  vs  x/b  for  fixed  time. 
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FIGURE   22.      Second-Order  V/ave    Form  Effects 
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VI.   CONCLUSIONS 

It  is  concluded  that  this  work  has  demonstrated  the 
feasibility  of  using  the  finite  element  method  to  solve  the 
linear  and  nonlinear  problems  of  the  type  presented.   Further, 
the  finite  element  method  is  shown  to  be  quite  versatile  in 
that  various  boundary  conditions  may  be  applied  to  the 
boundary  value  problems  studied,  with  only  minor  changes  in 
the  computational  scheme.   In  addition,  the  isoparametric 
elements  chosen  are  capable  of  representing  arbitrary  hull 
shapes . 

It  is  further  concluded  that  a  first-order  study  of  the 
effect  of  finite  depth  on  the  ship  motions  studied  demon- 
strates that  little  change  in  the  solutions  occur.   However, 
it  is  asserted  that  the  second-order  results  obtained 
demonstrate  that  nonlinear  (second-order)  effects  on  heaving 
motions  are  most  significant  in  very  shallow  water;  particu- 
larly for  motions  at  a  frequency  corresponding  very  closely 
to  the  natural  frequency  of  free  oscillations  of  the  floating 
body . 

The  present  state  of  digital  computer  capacity  appears  to 
preclude  the  extension  of  the  present  problem  to  three- 
dimensions.   However,  the  extension  of  this  study  to  other 
modes  of  motion  (sway  and  roll)  is  possible  as  demonstrated 
for  the  linear  sway  problem.   Further,  the  possibility  of 
studying  non-symmetric,  two-dimensional  hull  forms  is 
suggested. 
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An  adequate  finite  element  mesh  to  represent  the  region 
for  the  nonlinear  problems  studied  tends  to  tax  the  core 
capacity  of  the  computer  used.   Nevertheless,  meaningful 
results  were  obtained  through  careful  design  of  the  computa- 
tional scheme. 
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APPENDIX  A 
Complex  Algebra 

The  development  which  follows  determines  the  correct 
form  of  the  complex  numbers  w,  and  w?  in  the  following 
expression 


Re{z1eiat}Re{z2eiat}  =  Re^  +  w2e2iat}  (Al) 


z-.  and  Zp  are  complex  numbers  with  moduli  r-,  and  r?  and 
arguments  8-,  and  0p.   It  follows  using  appropriate  trigono- 
metric identies  that 


Re{r1e   1eiat}Re{r2e   2eiat}   = 


r  r?[cos61cose2  ^(l+cos2at)  +  sine,sin6p  ^(l-cos2at) 

■p  -p 

-  sin(01+62)^sin2at]  =-^?-2-   [cos(e1-02)    +    cos  (9-L+e2+2ct )  ]      (A2) 


Equation   A2   implies   that 

Z1Z2  Z1Z2  ,.-. 

w1   =  — p —  ,  w2    =   — 2 —      .  (A3) 

The  symbol  (~)  implies  conjugation. 

Equation  A3  shows  the  origin  of  the  factor  h   in  various 
places  in  the  text. 
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APPENDIX  B 
Isoparametric  Elements 
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FIGURE  Bl .   Element  Mapping  from  £,n  Plane  to  x,y  Plane 

The  essential  idea  in  the  isoparametric  element  lies 
in  the  following  equations  which  define  a  map  from  the  £,n 
plane  to  the  x,y  plane. 


*Te  e 
x  =  N  x    , 


y  =  Neye 


(Bl) 
(B2) 


re  . 


The  vector  N   is  a  1  x 12  row  vector  of  element  shape  functions 

e      e 
one  for  each  node  in  figure  Bl.   The  vectors  x  and  y   are 
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12  x 1  column  vectors  representing  the  respective  x  and  y 
coordinates  of  the  nodes  in  the  x,y  plane  (figure  Bl). 

The  shape  functions  are  conveniently  defined  using  the 
conventions  of  Zienkiewicz  [45].   Let 


«o  =  «i    >     % =  ™±        •  (B3) 


where  E,   and  n  are  the  coordinates  of  any  point  of  the  square 
element  in  the  £,r)  plane.   £.  and  n .  are  the  coordinates  of 
the  i   node.   The  shape  functions  are  given  by  the  following 
set  of  equations. 

Corner  Nodes: 


Ni  =  J2    (1  +  ^o)(1  +  %)["10  +  9(^  +  n2)]  '        (B2]) 


Edge  Nodes 


for  q  =  ±1,   n±  =   ±|  , 


N®  =  -^  (1  +  C0)(l  "  n2)(l  +  9n0)  ,  (B5) 


for  K±   =  ±|  ,    n±  =  ±1  , 


N|  =  ^  (l  +  n0)(l  -  C2)(l  +  9?0)  (B6) 
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The  element  is  defined  to  be  isoparametric  when  the  same 
shape  functions  are  used  to  achieve  the  element  mapping  as 
are  used  to  define  the  field  function  in  the  x,y  plane. 
There  are  certain  consistency  requirements  on  the  element 
shape  functions  as  discussed  by  Zienkiewicz  [45].   The  shape 
functions  chosen  satisfy  them  all. 

The  mapping  has  the  advantage  of  being  able  to  represent 
a  region  with  curvilinear  sides.   The  area  and  line  integrals 
given  in  part  III  are  calculated  by  the  expressions  given 
below. 

From  the  calculus,  the  following  relation  is  known  to  hold 


ft 

Ne 


~x 


r 


(B7) 


J  is  the  real  2x2  Jacobian  matrix  of  the  transformation 


,,e  e 
Nrx 


r?e  e 
N  x 


hi 


Me  e 

V 


(B8) 


It  follows,  provided  J  is  not  singular,  that 


Vx" 

Ne 

_~y 

.-1 


£ 


(B9) 


93 


«  \ 


Then,  for  example,  equation  86  in  part  III  may  be  written  at 

the  element  level  as 

,e 


T     T 

He  =  /   [N®     N®  ] 
z  Re   ~x     ~y 


N 
~x 


~y 


dRe 


(BIO) 


Using  equations  B8,  B9,  and 


dRe  =  dxdy  =  |J|  d£dn   , 


(Bll) 


equation  B9  may  be  written 


+1+1     T     T 

He  =   /   /   [N?     N®  ]  [J"1]7  [J"1] 
~     -1  -1    ~*>  ~n     z  ~ 


»t 


Jl  d£dn  .   (B12) 


Equation  B12  demonstrates  the  method  of  the  evaluation 
of  the  matrix  He  or  in  general  the  matrix  H.   The  integrations 
are  accomplished  using  six-point  Gauss  quadrature. 
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APPENDIX  C 
Computer  Program  Economization 


- 


FIGURE  CI.   Schematic  of  Coefficient  Matrix 


J  iCH 


The  economizations  in  the  computer  solution  of  the  linear 
system  of  equations  defined  by  equation  102  in  part  IV  are 
achieved  by  observing  certain  properties  of  the  matrix  K 
First  all  interior  nodes  of  a  given  mesh  contribute  only  real 
numbers  to  K    .   This  is  also  true  of  nodes  on  S.,,  Sp  and  S~ 
Only  the  nodes  on  Sj.  contribute  complex  entries  in  K    . 
Further  only  nodes  on  S~  and  Su    contribute  entries  which  are 
frequency  dependent. 

The  nodes  on  S„  and  S^   typically  comprise  no  more  than 
fifteen  percent  of  K^   .   Further,  the  nodes  on  S^  alone 
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comprise  nor  more  than  five  percent  of  the  total  rows  in  K    . 
Therefore  if  K    is  rearranged  with  equations  for  nodes  on 
Sjj  at  the  bottom,  equations  for  S^   nodes  directly  above  the 
Sn  equations,  and  all  other  rows  of  the  matrix  filling  the 
remaining  part  of  K    ,  the  resulting  matrix  would  appear  as 
schematically  represented  in  figure  CI.   The  region  designated 
with  a  Roman  numeral  I  (uncrosshatched)  indicates  the  portion 
of  the  matrix  completely  real  and  frequency  independent. 

Region  II  indicates  the  portion  of  the  matrix  which  has 
elements  that  are  frequency  dependent.   And,  Region  III 
indicates  the  portion  of  the  matrix  which  is  complex  and 
frequency  dependent.   If  the  matrix  K    is  formed  as  shown 
in  figure  CI,  eighty-five  percent  (approximately)  of  it  may 

be  eliminated  only  once  using  real  arithmetic.   This  is 

2 
possible  because  the  matrix  -a  QQ  does  not  have  to  be  added 

to  Region  II  until  after  Region  I  has  been  reduced  by  Gauss 

elimination,  and  the  same  fact  applies  to  the  matrix 


,4-CtJ 


(ia/cn)D  in  Region  III.   Then,  after  Region  I  is  eliminated, 

p 
the  matrix  -a  QQ  is  added  to  Region  II   (after  saving  the 

previous  results)  and  K    is  reduced  to  upper  triangular 
form  through  Region  II  using  real  arithmetic.   Now,  if  this 
result  is  saved  and  then  the  matrix  (ia/c-.)D  is  added  to 
Region  III  only  a  small  set  of  complex  equations  (usually 
about  25)  is  left  to  solve. 

After  the  final  elimination,  the  resulting  matrix  is  in 
complete  upper-triangular  form  and  ready  for  a  back- 
substitution  process.   The  right-hand  side  vector  is  complex 
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and  frequency  dependent  in  general.   However,  the  required 
information  to  perform  the  forward  substitution  process  on 
it  is  contained  in  the  final  form  of  K     (upper  triangular) . 
Therefore,  the  forward  substitution  on  ibar    may  be  done 
after  K    is  processed.   Then  a  back-substitution  of  the 
resulting  system  using  mixed-mode  computer  arithmetic  yields 
the  final  solution  vector  <jr  '  . 

A  close  examination  of  the  process  just  described  reveals 
that  for  a  whole  family  of  frequencies  only  a  small  portion 
of  K    has  to  be  reassembled  and  eliminated  again.   These 
matrix  properties  contribute  toward  a  sizeable  saving  in 
computer  time. 
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